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SUMMARY 


Computational procedures are presented for evaluating the 
sensitivity derivatives of the vibration frequencies and eigenmodes of 
framed structures. Both a displacement and a mixed formulation are used. 
The two key elements of the computational procedure are: a) Use of dynamic 
reduction techniques to substantially reduce the number of degrees of 
freedom; and b) Application of iterative techniques to improve the accuracy 
of the derivatives of the eigenmodes. The two reduction techniques 
considered in the study are the static condensation (Guyan's reduction) and 
a generalized dynamic reduction technique. 

Error norms are introduced to assess the accuracy of the eigenvalue 
and eigenvector derivatives obtained by the reduction techniques. The two 
iterative procedures used in improving the accuracy of the resulting 
eigenvector derivatives are based on the inverse power method and the 
preconditioned conjugate gradient technique. 

The effectiveness of the methods presented is demonstrated by three 
numerical examples. The structures considered are composed of beam elements 
and derivatives are taken with respect to cross-sectional areas and moments 
of inertia. An elastic line representation of a helicopter structure is 
also investigated and some of the concentrated masses are used as design 
variables. 
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I . INTRODUCTION 


1.1 Structural Sensitivity Analysis and Dynamic Reduction Methods 

There are an ever increasing number of fields where the ability to 
calculate the change in the response of a system due to a small change in a 
particular system parameter is making important contributions. A recent 
survey [1] reviews applications in such diverse disciplines as physiology, 
thermodynamics, physical chemistry and aerodynamics. Consequently the 
branch of science, known as sensitivity analysis, is the object of a great 
deal of research looking for better and more efficient ways to calculate 
sensitivity derivatives. 

The application of sensitivity analysis to structural response 
quantities, which has its roots in the development of gradient based 
automated structural optimization in the 1960's [2], is finding many new 
applications and gaining recognition as a powerful design tool [3]. 
Derivatives of static response, eigenvalues and eigenvectors, and transient 
response are being used in approximate analysis, for guidance in design 
modification and for improving analytical models. For extremely large and 
complex structures the calculation of the response and its derivatives can 
become prohibitively expensive unless some type of approximations are 
introduced. 

Today almost all dynamic response calculations of large structures 
are carried out using finite element methods. Due to the complex topology 
it is generally convenient to include much more detail in the model than is 
required to accurately predict the dynamic response [4]. This has resulted 
in the development of dynamic reduction methods that allow the number of 
degrees of freedom to be reduced without sacrificing accuracy. While quite 
a few methods have been developed for calculating eigenvalues and 
eigenvectors, no general and efficient method exists for calculating the 
derivatives from a reduced set of equations. The need to reduce the cost of 
calculating free vibration sensitivity derivatives is a strong motivation 
for applying dynamic reduction methods to this problem. 

1.2 Derivatives of Vibration Eigenvalues and Eigenvectors 

Free vibration eigenvalues and eigenvectors of large dynamic 
structures are frequently used in evaluating the response. Accordingly, 
their derivatives with respect to design variables provide useful 
information to the designer and analyst. For example, the approximate 
dynamic response of a structure that has been changed only slightly can be 
calculated quickly and easily. Also, derivatives of eigenvectors can be 
used as a guide to alter the response of a structure to minimize 
displacements at certain locations. Another application occurs in the 
development of analytical models where the derivatives of eigenvalues and 
eigenvectors can be invaluable for improving the model characterization to 
agree better with experimental results. 

The most straightforward method of calculating eigenvalue and 
eigenvector derivatives is by finite difference approximations. While this 
method is easy to implement, sensitivity to perturbation step sizes and the 
expense of analyzing the structure for at least one perturbation of each 
design variable are major drawbacks. Although recent progress has been made 
in the selection of step sizes [5], in many cases the method is 
computationally more expensive than other direct analytical methods. 
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The derivatives of eigenvalues can be calculated using a very simple 
analytical expression [6] that involves the eigenvalue being differentiated 
and its corresponding eigenvector. This is generally accepted as the 
preferred method. In contradistinction, the derivatives of eigenvectors can 
be calculated by several methods with the best choice being somewhat problem 
dependent. The modal method [6] calculates the derivatives as a linear 
superposition of the eigenvectors. While it is frequently possible to 
obtain accurate results using a truncated subset of eigenvectors, it is 
still computationally expensive and requires convergence checks. Other 
methods attempt to solve for the derivatives directly but must deal with the 
singularity of the governing differential equations. One procedure [6] 
performs algebraic manipulations that unfortunately destroy the banded form 
of the coefficient matrix; a major drawback when large matrices are 
involved. Another method [7] maintains the banded form of the coefficient 
matrix by calculating the derivative as the sum of a complementary and 
particular solution. In a majority of cases this method provides the best 
results; however, it does require the decomposition of a coefficient matrix 
for each distinct eigenvalue. A comparison of a number of methods for 
calculating vibration eigenvector derivatives is given in reference [8]. 

All the aforementioned methods for determining eigenvector 
derivatives provide direct f inite-number-of-step formulae. Reference [9] 
presents an iterative scheme that converges to the derivative of the 
eigenvector associated with the largest eigenvalue and suggests 
modifications that allow derivatives of additional eigenvectors to be 
calculated. The present study focuses on the efficient calculation of the 
derivatives of the eigenvectors associated with the smallest eigenvalues. 
However, the procedure is easily modified to converge to these vectors. 

1.3 Dynamic Reduction Methods 

There are practical engineering problems for which it is necessary 
to perform a dynamic analysis of very large structures by developing 
analytical models in terms of motions or forces at discrete points. A few 
examples include aircraft, rotorcraft, skyscrapers and nuclear power plants 
[4]. Before developing the model a decision must be made about the level of 
discretization desired. Models can be very coarse and lump many structural 
components together or extremely fine and include each structural component 
separately. The former case has the advantage of resulting in a much 
smaller number of degrees of freedom (hundreds as compared to thousands for 
the latter) and usually provides acceptable accuracy but masks much of the 
detailed information required about the behavior of the actual structural 
components. To get complete information on the dynamic stresses and 
displacements occurring in the structural components a very fine model is 
needed but this results in an extremely large number of degrees of freedom 
and high computational costs. This dilemma has resulted in the development 
of dynamic reduction methods which allow the size of the equations to be 
decreased without losing the fine detail of the model. 

The crux of dynamic reduction methods is the relationship between 
the small number of degrees of freedom defining the reduced model and the 
vector of unknowns of the original (full) model. This is accomplished by 
approximating the full vector of unknowns as a linear combination of a small 
number of global approximation (or basis) vectors, the amplitudes of which 
constitute the unknowns of the reduced equations. There have been a number 
of methods proposed for determining the basis vectors and a few of those 
will be described subsequently. 
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One of the first and most popular methods of dynamic reduction is 
known as static condensation or Guyan's reduction [10]. The procedure is to 
subdivide the vector of unknowns into two sets. One set contains the 
degrees of freedom considered most important that will be retained while the 
other contains those of lesser importance that will be eliminated. The two 
sets of degrees of freedom are usually referred to as free and constrained. 
The basis vectors are determined from the static relationship between the 
two sets and neglect the effects of the dynamic terms. This method is very 
easy to implement, however it requires some experience in choosing the free 
degrees of freedom, and in most cases only a moderate reduction is possible 
without sacrificing accuracy [4]. 

Dynamic condensation [11] is similar to static condensation but 
includes the dynamic terms in the relation between the free and constrained 
unknowns by using an iterative procedure. The method offers substantial 
improvements in accuracy, however because the equations must be reduced for 
each eigenvalue a great reduction in efficiency results (compared with 
static condensation). 

In both static and dynamic condensation the reduced unknowns are 
physical degrees of freedom. Another group of reduction methods uses a 
reduced set of unknowns which have no direct physical significance, thereby 
relieving the analyst of the task of choosing the free degrees of freedom. 
The first method of this type [12] proposed the use of Lanczos vectors, 
which are rich in the desired eigenvectors, for basis vectors. This has the 
advantage of automatic selection of basis vectors and can allow very 
substantial reductions in problem size. Subsequently, several variations on 
this original algorithm have been proposed such as subspace iteration [13] 
and generalized dynamic reduction [4]. These methods have become very 
popular in recent years due to their simplicity, accuracy and efficiency. 

1.4 Mixed Formulation Finite Elements 

In recent years there have been a number of applications in which 
the use of a mixed formulation in which the fundamental unknowns are both 
stresses and displacements has shown advantages over the more common 
displacement-based formulation. Studies in approximate structural 
reanalysis [14,15] have shown that the use of a mixed method provided higher 
accuracy for the predicted response, especially for stress quantities. 
Also, work in nonlinear structural analysis [16] has indicated that the use 
of reduction methods in conjunction with mixed finite element models can 
considerably improve the accuracy. 

When reduction methods are used in conjunction with the mixed 
formulation, explicit approximations are made for both the stresses and 
displacements. The stresses are allowed to be discontinuous across element 
boundaries and eliminated on the element level. It has been shown [16] that 
this improves the performance of the method and also avoids the 
computational expense that would be associated with generating the basis 
from the complete set of equations directly. 

1.5 Review of Previous Pertinent Work 

While several examples of the current work being performed in the 
fields of free vibration structural sensitivity analysis and dynamic 
reduction methods have been discussed previously, there are surprisingly few 
cases in which the two are considered jointly. A recent paper on the 
calculation of eigenvalue sensitivity derivatives [17] proposes two new 
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solution techniques, one direct and the other iterative, and also addresses 
the case of repeated eigenvalues. However, no reduction in problem size or 
significant improvements in efficiency are suggested. Another study [18] 
suggests improvements to increase the rate of convergence of the modal 
method; however, it was found in [8] that even with these improvements the 
method is computationally expensive. 

A paper by Robinson [19] describes the implementation of a 
systematic finite element modification technique in the program EAL for the 
purpose of improving dynamic models to agree better with experimental 
results. In this implementation the eigenvector derivatives are calculated 
by the modal method and he suggests the use of a Rayleigh-Ri tz type 
reduction method to reduce the computational expense. In the described 
implementation, the modification of the structure is performed in a reduced 
space defined by the eigenvectors of the original structure. Large 
reductions in computer execution times are reported. 

Another important application of free vibration sensitivity 
derivatives is in automated structural optimization. Kim and Anderson [20] 
apply generalized dynamic reduction to automated redesign for frequency 
changes. The optimization is performed in the reduced subspace and good 
results are reported; however, changes in eigenvectors are not considered. 
Also, the optimization equations are formulated in terms of perturbation 
equations that do not require the calculation of gradients, but this study 
does give an indication of the ability of the subspaces to accurately 
represent the full dynamic space. 

The calculation of approximate eigenvalue and eigenvector 
derivatives for a modified system is addressed in a recent paper by Wang 
[21]. For the special case in which the stiffness and mass terms are 
homogeneous functions of the design variable he develops expressions using 
cross mode energies to calculate the approximate derivatives. The results 
that are reported show a mixed degree of success. 

1.6 Objective and Scope of the Present Study 

The objective of the present study is to develop an efficient 
computational strategy for calculating the sensitivity derivatives of 
vibration frequencies and eigenmodes of framed structures. The key elements 
of the strategy are: 

1. Use of dynamic reduction methods to substantially reduce 
the number of degrees of freedom. 

2. The application of a mixed method in conjunction with 
reduction methods to the free vibration derivative 
calculations . 

3. Use of iterative techniques for improving the accuracy of 
the derivatives obtained by using the reduction technique. 

The scope of this study includes: 

1. Large frame-type structures modelled by finite elements. 

2. Sensitivity derivatives of vibration eigenvalues and 
eigenvectors with respect to cross-sectional design 
variables and concentrated masses. 
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II. SENSITIVITY CALCULATIONS VIA REDUCTION METHODS 


2.1 Introduction 

Reduction methods is the name usually given to a class of hybrid 
two-step computational procedures that attempt to reduce the dimension of a 
given problem and permit the solution to be calculated more economically. 
The first step is the spatial discretization of the system into a number of 
finite domains within which interpolation functions are used to approximate 
the unknown fields in terms of unknown nodal parameters. These nodal 
parameters span what will henceforth be referred to as the full space of the 
system. The second step is to impose relationships between the degrees of 
freedom of the full space to define a new coordinate system, usually of much 
smaller dimension, called the reduced subspace. The transformation of the 
governing equations from the full space to the reduced subspace is 
accomplished by using a Rayleigh-Ri tz (or Bubnov-Galerkin) procedure. 
Throughout this report the words full and reduced are used as a prefix to 
indicate the space in which the response quantities are defined. For 
example, reduced equations refer to equations defined in the reduced 
subspace. 

In this section the derivation and solution of the reduced equations 
for the derivatives of vibration eigenvalues and eigenvectors is discussed. 
Since this is accomplished by a transformation of coordinates the full 
equations are developed first and are then transformed into the reduced 
subspace. 

2.2 Governing Equations 

The free vibrations of structures modelled by displacement-based 
finite elements are governed by the following system of equations: 

( [K] - X.[M] ) {X.} = 0 (la) 

where 

[K] = global stiffness matrix 

[M] = global mass matrix (consistent or lumped) 

X^ = ith eigenvalue (square of the frequency of vibration) 

{X^} = ith eigenmode (vibration mode shape) 

The eigenvectors satisfy the orthonormality relation 

{X.} T [M] {X.} = (lb) 

where 

5. . = Kronecker delta =1 if i=j 
ij J 

= 0 if i*j 

Before assembling the stiffness and mass matrices the specific 
properties of each element must be known. Changing any one of these 
properties would change the resulting stiffness and mass matrices and thus 
the eigenvalues and eigenvectors. To quantify the effect of this change the 
equations (1) are differentiated with respect to the property, called a 
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design variable [6]. It is often more convenient to determine the effect of 
the simultaneous change in the properties of several elements and this is 
accomplished using design variable linking. The properties of each 
individual element (dependent design variables v ^ ) are related to a much 

smaller number of parameters (independent design variables V^) by a linking 

relationship. For example, linear design variable linking is expressed by 
the relation 


v. = E c., V. 

■j k=i Jk k 


where 


n = the number of independent design variables 
c^ = the linking coefficients 


( 2 ) 


Differentiating equations (1) with respect to an independent design variable 

V. results in 
k 


and 


( [K] - 
2 {X.} T 


9{X } 

VM] > -9V 1 - 


3{X.} 

+ 



(3a) 

(3b) 


The unknowns in equations (3) are the derivatives of the eigenvalue and 
3X 3{X } 

eigenvector, -rrj- and —zr . — . It is noted that when using design variable 

9V k 9V k 

linking the derivatives of the stiffness and mass matrices must be evaluated 
using the chain rule of partial differentiation. For the linear design 
variable linking in equations (2), this results in 


9[K] 

m 

= E 

9[K] 

9v . 

_I _ 

m 

l 

9[K] 

c . . 

9V k 

j=i 

9v . 

3 

1 

> 

ro 

j=i 

9v . 

3 

jk 

9[M] 

m 

= E 

9 [ M ] 

*1 

m 

3[M] 

c . , 

9V k 

j=l 

9v. 

J 

1 

> 

CO 

j=l 

9v . 

3 

jk 


(4a) 

(4b) 


where m is the number of dependent design variables. 

An examination of equations (3a) indicates that a direct solution is 
not possible for two reasons: 


1. The right hand side vector contains the derivative of the 
eigenvalue which is, as yet, unknown. 

2. The matrix on the left-hand side is singular. 
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Before solving equations (3) it is important to determine the multiplicity 
of the eigenvalue X^ . Since the solution for nonrepeated eigenvalues is 

simpler it is presented first and then generalized to the repeated case. 
2.2.1 Nonrepeated Eigenvalues 


The derivative of the eigenvalue is calculated through premultiplication of 
equation (3a) by the transpose of the eigenvector and noting equations (1), 
resulting in the following expression: 


3X. 

* 



3[K] 

^iT 


3[M] 


X i 3V, 


) (X^ 


(5) 


Nelson [7] solved the problem of a singular coefficient matrix by expressing 
the derivative of the eigenvector as the sum of a particular and 
complementary solution. 


3{X } 

-svr " ">i> + c 'V 

k 


( 6 ) 


The complementary solution {X^} is the eigenvector and the particular 

solution (Q. } is obtained by constraining one degree of freedom of the 

solution vector in equation (3a) to zero (e.g. the one corresponding to the 

largest component of the eigenvector) to eliminate the singularity in the 

coefficient matrix. The unknown constant C is then determined by 

substituting the assumed form of the solution (6) into equation (3b). This 
leads to 


C = - 


{ V 


[M] {Q.} 


1 T 

i (V 


3 [ M ] 




(7) 


2.2.2 Repeated Eigenvalues 

When eigenvalues of multiplicity greater than one are present the 
solution of equations (3) becomes much more complex due to the nonuniqueness 
of the associated eigenvectors and the decreased rank of the coefficient 
matrix. For simplicity an eigenvalue of multiplicity two is assumed in the 
subsequent discussion however, the method is easily generalized to cases of 
higher multiplicity (see for example Reference [17]). The free vibration 
equations for the repeated eigenvalue can be written 

( [K] - A.[M]) {X.} = 0 (8a) 

([K] - X. +1 [M]) {X. +1 } = 0 (8b) 

where 
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The eigenvectors satisfy the orthonormality relation (lb). Since any linear 
combination of eigenvectors is also an eigenvector the free vibration 
equation can also be expressed as 


UK] - X.[M]) {XC} = 0 (9) 

where 

{X i } = tX i X i + 1 ] {a i } < 10 > 

The matrix [ X ± X i+1 ] stores the eigenvectors columnwise and {Xp also 

satisfies the orthonormality relation (lb). Differentiating (9) with 
respect to an independent design variable results in 

3{Xn 3X 3[M] 3[K] 

UK] - V M » Tv^ - (w^ |M l + x i^ - T^J < n > 


Premultiplication of (11) by [X^ and noting equations (8) and (10) 

leads to the following: 


[X. X 

1 l 


i+1 



9[M] 

X i-3V^J I X i X i + ll < a i> 


ax. 



( 12 ) 


Equation (12) is an eigenvalue problem of dimension two and has the two 
3X. 3X i + l 

solutions {a^} and , {a^} which indicates that the repeated 

k k 

eigenvalue may follow two separate paths with the continuous variation of 

and also leads to two separate right hand side vectors for equation (11) 
(and thus two eigenvector derivatives). Differentiation of the 
orthonormality relations for the vectors {X!} and {X( + ^} (obtained from {a.} 

and {a^ + ^} respectively) results in the following equations: 


t 3{XU , T 3[M] 

ixp [H] -5v- - - 2 (xp -jy- (xp 


{Xi +1 } A [M] 


3 < x i + l> 

3V,_ 


1 T 3[M] 

2 {x i+i } “avT 


< x i + i> 


(13a) 

(13b) 


mU 1 t 3[M] t 3{X' } 

-W- l M HXC +1 } + } + (Xt j [M] = 0 

k k k 


(13c) 


As in the previous subsection, solutions to equation (11) are assumed of the 
form 


3 (XU 

ivr = «i> + i x i x i + ii 

k 


(14a) 
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3 ' x i + l> 


(14b) 


av, 


■ «ui> * l x i x ui' <“> 


In equations (14), {Q^} and {Qj, + ^} are particular solutions to (11) obtained 

by constraining two degrees of freedom (those associated with the largest 
degree of freedom of each eigenvector) in the solution vector and the 
eigenvectors are complementary solutions. Substitution of the solutions 
(14) into equations (13) leads to the following expressions for the unknown 
constants : 


, t 3[M] t 

i - -5 < X P SvT {x i> - f x i} IMKQi) 

,3[M] 


c, = 


i ' x U> T -WT {x i.i ! - f x i + i) T l H lW it i) 


av,. 


,3[M] 


' 2*1 * - (Oi> T [«H x i +1 > - ( x -! T [«IW i+1 ) - ( x i) T ^r( x i + l! 

k 


(15a) 

(15b) 

(15c) 


It is noted from equation (15c) that the nonuniqueness of the eigenvectors 
extends to their derivatives. 

Equations (3) or (11) can be solved by several alternate methods 
(see subsection 1.2) but the solution presented herein has the advantage of 
providing essentially exact results without the need for convergence checks. 
Its major disadvantage is that the coefficient matrix must be decomposed for 
each distinct eigenvalue. However, it is shown subsequently that reduction 
methods can substantially reduce the size of this matrix and the cost of 
decomposition. 

2.3 Dynamic Reduction Methods 

Dynamic reduction methods approximate a vector in the full space by 
a linear combination of a small number of linearly independent basis vectors 
that define a reduced subspace. The accuracy of the method depends on the 
degree to which the full vector lies within the reduced subspace and, as 
such, the determination of the basis vectors is a crucial step that depends 
strongly on the nature of the full vectors. Since the current work focuses 
on the derivatives of the eigenvectors, a logical starting point is methods 
developed for approximating the eigenvectors. The approximation is 
expressed by the following equation: 

{X.} - {X.} = [T] {*.} (16) 

where 

{X^} = the full eigenvector 

{X.} = the approximate full eigenvector 

in = matrix with basis vectors as columns 

{^} = reduced eigenvector 
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The reduced eigenvector consists of generalized coordinates that are 

amplitudes of the basis vectors. 

Using equation (16) the free vibration equations (1) are transformed 
into the reduced subspace resulting in the following reduced equations; 


( IK] - X. [M] ) {<!/.} = 0 (17a) 

and 

{*.} T [M] {^} = Sjj (17b) 

where 

T 

[K] = [T] [K] [T] = reduced stiffness matrix 

T 

[M] = [T] [M] [T] = reduced mass matrix 

Xj = approximate eigenvalue 
{ ) = reduced eigenvector 

Since equations (17) represent a much smaller system of equations than (1), 
their solution requires substantially less effort and provides 
approximations to the lower eigenvalues and, through the use of equation 
(16), the eigenvectors. Several methods have been developed for performing 
reduced eigenvalue analysis (see for example [4] and [11]) and two of these 
are discussed in the following subsections. 

2.3.1 Static Condensation 


One of the first and most popular methods of determining the basis 
vectors used in equations (16) and (17) is known as static condensation or 
Guyan's reductions [10]. The procedure begins with the choice of 
constrained degrees of freedom (a subset of the degrees of freedom of the 
full space that are assumed massless). Then the equations (1) are 
partitioned in the form 


/ 

. K ff K fc . 


• M ££ M fc ■ 

\ 

V 

. K c£ K CC . 

- X. 
1 

. h c£ H CC . 

/ 



(18) 


where the superscript f denotes the free (retained) degrees of freedom and 
the superscript c denotes the constrained degrees of freedom. Expanding the 

equation in the lower partition of (18) and solving for {X?} results in 

{X?} = - ( [K CC ] - X. [M CC ] T 1 ( [K cf ] - X. [M cf ] ) {X*} (19) 

c c c f 

If the mass matrices [M ] and [M ] are neglected, equation (19) reduces to 
the following static relation; 

{X^} = - [K CC ] _1 [K cf ] {X^} (20) 
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Equations (20) are used to eliminate the constrained degrees of freedom. 
The matrix of basis vectors of the static condensation procedure can then be 
expressed by the equation 


m = 


[i] 


[K CC ] 1 [K cf ] 


( 21 ) 


cc c f 

The approximation results from neglecting the mass matrices [M ] and [M ]. 
The accuracy of the approximation depends on the selection of free and 
constrained degrees of freedom. There is obviously a certain amount of 
talent or experience required for the proper selection of the free and 
constrained degrees of freedom. Also, looking at equation (19) it is 
apparent that the effect of the neglected masses increases with increasing 
eigenvalue and therefore, the accuracy of the method decreases with 
increasing eigenvalue. 

With static condensation it is generally possible to perform a 
moderate amount of reduction and still obtain good results for the lowest 
few eigenvalues and eigenvectors. However, the desire to obtain an even 
higher degree of reduction and also to eliminate the need for choosing 
massless degrees of freedom has led to the development of another group of 
methods that generate the basis vectors automatically. These are discussed 
in the succeeding subsection. 


2.3.2 Generalized Dynamic Reduction 


Since the main purpose of the dynamic reduction methods is the 
accurate calculation of eigenvectors, procedures have been developed that 
start with arbitrarily chosen vectors and then attempt to enrich the content 
therein of a selected group of eigenvectors (see for example [12] and [13]). 
The enriched vectors, called Lanczos vectors, are then used as basis 
vectors. In generalized dynamic reduction [4], the Lanczos vectors are 
generated by performing a required number of iterations through the 
recursion relation 

( [K] - v[M] ) [ Y] (r+1) = [ M ] [Y] (r) (22a) 

where 


m 


(r) 


= [{Yj} {Y 2 } ... {Y n }] 


and 


. (r+1) = (r+1) ( y ,(r+l) 

1 1 i i l 1 i J 


(22b) 

(22c) 


In equations (22), the constant a^ r+ ^ is used to normalize {Yj}^ r+ ^1 


As 


mentioned previously, the initial iteration vectors [Y] are chosen 
arbitrarily. Experience has shown that the use of a pseudorandom number 
generator (e.g. generated by the multiplicative congruential method, see 

Reference [22]) in initializing [Y]^^ helps insure that all eigenvectors 
are represented and increases the chances that a mode will not be missed. 
In addition, if a structure contains large concentrated masses it has been 
found advantageous to include unit vectors to excite these masses. Since 
eigenvectors other than the ones associated with the smallest eigenvalues 
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are sometimes of interest, a nonzero value for the shift parameter v allows 
equation (22) to converge to the eigenvectors desired. 

In addition to being rich in the eigenvectors of interest the basis 
vectors must also be linearly independent. The final iteration vectors from 
equation (22) are orthonormalized by performing the Gram-Schmidt 
orthogonalization procedure. Often this must be performed several times to 
achieve the high degree of linear independence required to insure the 
solvability of equation (17). The orthogonalization also helps to amplify 
the content of eigenvectors slightly further away form the shift value \>, in 
the basis vectors. 

Generalized dynamic reduction usually allows a greater reduction in 
the number of degrees of freedom than static condensation and at the same 
time provides greater accuracy. The number of modes that can be obtained 
accurately from the reduced equations is usually less than one-half the 
number of basis vectors. With static condensation this limits the number of 
modes that can be calculated, but generalized dynamic reduction used with 
multiple shift values v and multiple reductions can be used in calculating 
many more eigenvectors than the number of basis vectors used. 

2.4 Reduction of the Governing Finite Element Equations 

The preceding section has dealt with dynamic reduction methods for 
calculating eigenvalues and eigenvectors. The next step is to apply these 
to the calculation of the required derivatives. Differentiating equations 
(16) and (17) with respect to an independent design variable results in 


and 


where 


and 


3{x } 3{x. ) 9{* } am 

av k ~ av k = av k + 1v^ **i* 

3{* } aX a[M] a[K] 

( [K] - A. [M] ) -gy— ( l m 1 + \ “avjj _ “av^ ) 

3[K] m 3[K] m 3 [ T] 3[ r] T 

“IvT = I r l -avT ^ + ™ W ^vT + ^vT [K] [r] 

k k k k 

3[M] 3[H] am 31 r] T 

Tvr - rn T -5T in * m i»i ir * it i«i m 

k k k k 


(23) 


(24) 


The matrix 3[r]/3V k in equations (23) and (24) represents the changes in the 

basis vectors due to changes in the design variables. It is not desirable 
to calculate this matrix since it would be at least as computationally 
expensive as calculating the derivatives of the eigenvectors from the full 
equations. The purpose of this procedure is to reduce the expense and the 
contribution of the derivatives of the basis vectors must be neglected to 
retain the desired efficiency. With this assumption the reduced derivative 
of the mass and stiffness matrices from equation (24) become 
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(25a) 


and 


a[R] 

9[M] 

~W, 


= m 


= [r] 


a[K] 

3[M] 

“9VT 


m 


m 


k k 

and the derivative of the eigenvector is approximated as 


(25b) 


3{X.} 9{X.} 


m 


3{^} 


(26) 


It is noted that neglecting the derivative of the basis vectors is 
equivalent to assuming the derivative of the eigenvector lies in the 
subspace defined by the basis vectors. 

Since the form of equation (24) is identical to that of equation (3) 
only of much smaller order, it can be solved by the same method. For 
eigenvalues of multiplicity one the approximate derivative of the eigenvalue 
is calculated from 


3X. - 3[K] 9[M] 

av^ = **1* ( \ “9V^ ) { ^i } 


(27) 


The reduced derivative of the eigenvector is assumed to be of the form 


St'J'- } 

~mr - < 9 i> * c «»i> 

k 


(28) 


where the reduced eigenvector {\J^} is a complementary solution and {9^} is a 

particular solution to equation (24) obtained by constraining one component 
of the solution vector. In the reduced equations the degree of freedom 
associated with the smallest term on the diagonal of the coefficient matrix 
is constrained to zero. The constant c is then calculated from 


c = — 


IM] {9.) - \ {*.} T 


3[M] 

~9VT 




(29) 


For repeated eigenvalues a procedure analogous to that used above, but using 
the method described in subsection 2.2.2 can be used. 

The foregoing procedure for calculating approximate derivatives of 
eigenvalues and eigenvectors with respect to design variables greatly 
reduces the computational effort compared to the solution of the full 
equations. The performance of these methods on several finite element 
models is discussed in section five in conjunction with the numerical 
studies. 
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III. APPLICATION OF A MIXED FINITE ELEMENT MODEL 
TO SENSITIVITY CALCULATIONS 


3.1 Introduction 

In the preceeding section reduced equations are developed for 
calculating the derivatives of free vibration eigenvalues and eigenvectors 
based on a displacement formulation. The displacement formulation is the 
most commonly used; however, in some applications a mixed formulation, which 
takes both stresses and displacements as unknowns, has several advantages 
(see subsection 1.4 for references). Studies in nonlinear analysis have 
shown that the use of a mixed model in conjunction with reduction methods 
can provide improved accuracy over reduction methods used with a 
displacement formulation [23], Also, in several types of structures (e.g. 
trusses, plates, membranes) approximate reanalysis is improved by taking the 
inverse of the sizing quantities as design variables [24]. In the mixed 
formulation this is a natural choice. 

3.2 Element Formulation 


The mixed finite element formulation is based on the Hellinger- 
Reissner variational principle [25]. The following functional, which 
represents a modified version of Hamilton's principle, is used in the 
development of the governing equations: 


" - j ' 2 ( “ - < v - V ) dt 

l l 

where 

V = f {e} T {t} d S2 

“ Q 

U = j \ {t} T [a] {t} da 
J Q 

and 

v = h f P [u} T {U> da 

J a 


(30) 


The integrations are performed over the volume of the structure Q. The 
infinitesimal strains {s} correspond to the displacements {u} of the 
structure from its undeformed state and the stresses {x} result from these 
strains. The matrix [a] contains the material compliance coefficients. It 
is noted that no external forces are included in the functional since the 
present study is considering only free vibrations. If considered, body 
forces, surface tractions and point loads would be added to the term W. 

The functional in equation (30) is discretized by subdividing the 
structure into finite elements and approximating the fundamental unknowns 
within each element using Lagrangian interpolation functions. The 

displacements and stresses of the j** 1 element are given by the relations 
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(31) 


{u j } = [N j ] {U j } 

{-r 5 } * [N j ] {H j } (32) 

where {U} are the nodal displacements and {H} are stress parameters. It was 
shown in [23] that better accuracy is obtained from mixed elements when the 
interpolation functions for the stresses are one degree lower than those 

used for approximating the displacements. Accordingly, [N] represents a 
matrix of lower degree shape functions than [ N ] . Also, the stress field is 
allowed to be discontinuous across element boundaries so that the stress 
parameters can be eliminated on the elemental level. 

Applying the aforementioned finite element approximation results in 
the following discretized functional: 


n = J tz ( w - (V - u c ) ) dt 

t l 

where 

V = l [ , {U} T [ B ] T [N] {H} dQ 

j=i v 

0=1 j [ . {H} T [N] T [a] [N] (H) dQ 
j=i Jq 3 

and 

V = S \ f . P {U) T [N] T [N] {U} dQ 

j=i J & 


(33) 


The summations are over the number of elements L, and [B] is the discrete 
strain-displacement matrix. 

Varying the displacement and stress parameters independently and 
simultaneously results in discrete governing equations which consist of both 
constitutive relations and equilibrium equations. For the j th element these 
equations are cast in the following form: 



(34) 


The submatrices in these equations are: 

[ F ] ^ = [ . [N] T [a] [N] dQ = elemental flexibility matrix 

J Q 3 

[S] 3 = I . [N]^ [B] dQ = elemental generalized stiffness matrix 
J O 3 

[M] 3 = [ . p [N] T [N] dQ = elemental consistent mass matrix 

J Q 3 

For free vibrations the nodal displacements (U) are of the form 
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{X(G,X)} 


(35) 


(U(S,t)} = e lAt 
and equations (34) become 


'0 0 1 J ( H. 'I 3 

1 = 0 (36) 

.0 M J l X. J 

± l 

Since the stress field is allowed to be discontinuous across interelement 

boundaries, the upper partition of (36) can be solved for {fP} at the 
element level. This allows the elimination of the stress parameters from 
the lower partition. The resulting elemental equations can then be 
assembled in the form 



( [K] - X. [M] ) {X.} = 0 (37) 

where 

L . 

[K] = Z [S J ] [F^] - [S' 1 ] = global stiffness matrix 

j=l 

and 

L 

[M] = Z [M^ ] = global consistent mass matrix 
3=1 

The above system of equations is of the same order as a corresponding 
displacement model. 

3.3 Derivatives of Eigenvalues and Eigenvectors 

Equation (37) is similar in form to the governing equations of the 
displacement formulation (equation(l)) and the sensitivity derivatives of 
the eigenvalues and eigenvectors can be calculated as outlined in equations 
(3) — (7) . The only difference is the derivative of stiffness matrix which in 
this case is 


3[K] 

~"k 



[s¥ 


aUF 3 ] 1 ) 
3VT 


[S J ] 


(38) 


This approach will work, but has no advantage over the displacement 
formulation since it requires the differentiation of the inverse of the 
flexibility matrix (essentially a stiffness matrix). In cases where it is 
advantageous to choose the inverse of the sizing variables as design 
variables (discussed at the beginning of this section) it would be more 
efficient to differentiate the flexibility matrix and thus an alternate 
formulation of the derivative equation is required. Differentiation of 
equation (36) with respect to an independent design variable results in 

the following: 


17 


-F S 
oT 


0 


- X. 
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0 0 




8V k 
3X. 

I 3V^ j 


ax. 

i 

3V, 



j 
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j 

r sf i 




0 0 


0 0 


0 



H. 


+ X. 



3v k 


< 

l 

0 M 

i 

o - M 

3V. 

L k J 


0 0 



l 

L / 


( 39 ) 


The solution of the upper partition of (39) for the derivative of the stress 
parameters, again on the element level, and the solution of the upper 
partition of (36) for the stress parameters can be used to eliminate all 
stress quantities from the lower partition of (39). The assembled global 
equations then take the form 


~ 3{X.} 

<[K] - X. [M] ) 


3X - 3 [ M ] 3 [ K ] 

( 3VT lM1 + X i“3VT + "IvT ) < X i> 


(40) 


where 


3[K] 


3V. 


= Z [ ] T [ ] 1 
k j=l 


3[F J ] . 


3V, 


[F^J 1 [S J ] 


The right hand side vector contains the derivative of the flexibility matrix 
as desired. Equation (40) is in the general form of equation (3) and can be 
solved by any of the several methods developed for calculating eigenvalue 
and eigenvector derivatives. 


3.4 Reduction Methods 


The basic idea behind the use of reduction methods in conjunction 
with mixed models is to define a small number of parameters related to both 
the displacement and stress degrees of freedom of the full model and reduce 
the equations to these coordinates. This is accomplished by defining a 
basis vector representing each reduced coordinate. The stresses and 
displacements are then approximated as a linear combination of the linearly 
independent basis vectors. In partitioned matrix format this approximation 
is 


H 

U 




(41) 
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where [ T R ] and [T^] are matrices with the stress and displacement basis 

vectors as columns, respectively, and {<f>} are the reduced unknowns. A major 
advantage of the mixed model is the direct approximation of the stresses. 

The basis vectors are determined in a two-step procedure. The 
displacement basis vectors [T^] are determined first by one of the methods 

used in section two. In static condensation, some of the displacement 
degrees of freedom are assumed massless and eliminated (see section 2.3.1). 
For the mixed model this leads to the following basis vectors: 


i r x ] = 


where 


[K] 


[I] 

— [ K CC ] ~ 1 [ K cf ] J 


K ff K fc ' 
K cf K CC . 


(42) 


is a partitioned form of the stiffness matrix defined in equation (37). The 
superscript f denotes degrees of freedom with mass and the superscript c 
those without. 

Generalized dynamic reduction can also be used to determine the 
displacement basis vectors (see section 2.3.2). Initial iteration vectors 

are chosen arbitrarily and placed as columns in the matrix [Z]^\ A 
required number of iterations is performed through the recursion relations 


where 


and 


( [K] - v [M] ) [Z] (r+1) = [M] [Z] (r) 
[z] (r) = l {Z 1 J {Z 2 } ... {Z n }] (r) 


{Z.} (r+1) = 3< r+1) {Z.} (r+1) 
1 1 1 1 


(43a) 

(43b) 

(43c) 


where the matrices [K] and [M] are defined in equations (34) and (37) and 

^(r+l) usec [ t 0 normalize {Z. }^ r+ ^ . The final iterates are 

1 1 i 

orthonormalized to form a linearly independent set. 

The second step is the determination of the stress basis vectors. 

As mentioned previously, for efficiency the matrices [F] and [S] in equation 

(34) are only calculated at the element level and never assembled globally 

except for an element by element contribution to the stiffness matrix as 

defined in equation (37), For this reason the stress basis vectors are also 

calculated only at the element level. The stress and displacement degrees 

of freedom in the jth element are related by the equation 

[F j ] {H j } = [S j ] {U j } (44) 

Using the approximation given in equation (41) the above relation can be 
solved for the stress basis vectors to yield 
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(45) 


[rj] = [F j ] _ 1 [s j ][rj] 

where [T^] is a matrix containing the rows of [1^] corresponding to the 

displacement degrees of freedom of the jth element. 

Equation (41) is substituted into the discretized functional 
(equation(33) ) and the resulting Rayleigh-Ritz analysis leads to the reduced 
equations 


( [K] - \ [M]) {*.} = 0 


The reduced matrices are given by 
L 


[K] = z [rj] T [F j ][rj] + [rj] T [s j ][rj] + [rj] T [s j ] T [rj] 

j=i 

= reduced stiffness matrix 

T ~ 

[M] = [ r x ] [ m ] [ r x ] = reduced mass matrix 


(46) 


(47) 

(48) 


where [F~*], [S^] and [M] are defined in equations (34) and (37) and it is 

noted that [K] is formed from element contributions. Also, for free 
vibration, the reduced unknowns {$} are assumed to be of the form 


{<f>(r,t)} = e iAt {^(T,X)} (49) 

The reduced equation for the eigenvalues and eigenvectors (46) looks very 
similar to the reduced equation (17) from the displacement formulation with 
the exception of the reduced stiffness matrix which has a very different 
form. 

A reduced equation for the derivative of the eigenvalues and 
eigenvectors can be developed by differentiating the reduced equation (46) 
with respect to an independent design variable V^> resulting in the 

following: 


where 



(50) 


The basis vectors in equation (42) are assumed to be invariant with respect 
to the design variable which, as discussed previously, is necessary to 
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obtain the desired efficiency. The advantage of this formulation over the 
displacement formulation results from the use of the derivatives of the 
flexibility coefficients on the right hand side (as opposed to the 
derivatives of the stiffness coefficients in (24)). Equations (46) and (50) 
can be solved by any of the methods described previously, and the 
eigenvector derivatives recovered from the relation 



(51) 


Numerical results obtained by using the mixed formulation are discussed in 
section five and compared with those from the displacement formulation. 
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IV. ACCURACY AND ITERATIVE REFINEMENT 


4.1 Introduction 

Methods for calculating the derivatives of eigenvalues and 
eigenvectors from reduced equations are developed in the preceding two 
sections. The usefulness of the proposed procedure requires the ability to 
assess the accuracy of the eigenvector derivatives and a method for 
improving the accuracy of the derivatives that do not meet the required 
tolerance. In order to accomplish these two objectives it is necessary to 
return to the full equations. Two methods are presented that use iterative 
techniques involving a shifted stiffness matrix to improve the eigenvector 
derivatives. As a by-product of the performance of generalized dynamic 
reduction this shifted stiffness matrix is decomposed and as such, the 
implementation of these methods adds minimal computational expense. 

4.2 Solution Errors 

A key element of approximate solution methods is the ability to 
evaluate the accuracy of the approximate solution without knowing the exact 
solution. In the present study this is accomplished by calculating a 
normalized residual error from the full equations. For example, the 
accuracy of the approximate eigenvalues and eigenvectors is evaluated using 
the following error norm [13]: 


II <[K] - X[M]){X) || 

s. = (52) 

» [KHX) || 2 

The numerator is the Euclidean norm of the residual vector from equation 
(la) and the denominator is a normalization factor (the norm of the internal 
force vector). Although this does not provide a pointwise error measure, it 
does gives a good indication of the overall accuracy. Similarly, the 
accuracy of the derivatives of the eigenvalues and eigenvectors can be 
evaluated using the following norm: 


( [K] - X[M] ) 


9 (X) raX 


£ 2 = 


The numerator is the 
denominator is the 
tolerances, s. or s„ 


3V 


(_[«! 

V av 


x3[M] _ 3[K1] { x } 
av av ' 


(53) 


1 ( 9X [M] 

+ X 3[M] _ 3[K]' 

1 {X} 

1 


11 V av 

av av ; 

f 

*2 


norm of 

the residual 

from 

equation (3) 

and the 

norm of 

the right hand 

side. 

If either of 

the error 

is not met 

then the corresponding approximate 

solution 


must be improved. 
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4.3 Improvement of the Eigensolution 

The error norm gives an indication of the accuracy of an 

eigenvalue and its corresponding eigenvector. If the required tolerance is 
not met then two possibilities are considered for improving the accuracy. 
The first is to improve the characterization of the reduced basis and solve 
for a new eigensolution (this results in increased accuracy for all the 
eigenmodes) or alternately, the particular eigenmodes not meeting the 
tolerance can be individually improved. When using static condensation, the 
former approach essentially involves starting over but with the benefit of 
additional knowledge. Improving the basis requires including the mass 
associated with more and/or different degrees of freedom. The approximate 
eigenvectors from the previous attempt are useful in this choice. In 
contradistinction, generalized dynamic reduction is much more amenable to 
improvement. The characterization of the reduced basis is augmented by 
performing more iterations on the existing basis vectors (or on the 
approximate eigenvectors) and also with additional vectors. In either case 
none of the information in the current basis is discarded. In addition, if 
only a few of the modes are found to be inaccurate, they alone are improved 
by cycling through the relations 

( [K] - v[M] ) {X} (k+1) = [M] {X} (k) (54) 

X (k+1) = {X} (k+1) [K]{X) (k+1) / {X} (k+1) [M]{X} (k+1) (55) 

which represent the inverse power method. Since the coefficient matrix in 
(54) has been decomposed previously, these procedures add little 

computational expense. The previous discussion brings forth some of the 
advantages of generalized dynamic reduction over static condensation. Due 
to its adaptability to the type of analysis being performed only generalized 
dynamic reduction is considered in the remainder of this section. 

Before turning to methods for improving the accuracy of the 
derivatives of the eigenvectors it is useful to identify the major time 

consuming operations in the calculation of the eigensolution. The largest 

percentage results from the decomposition of the shifted stiffness matrix 

used in the generation of the basis vectors. Two other items that account 
for a large percentage of the computer time are the back substitutions 
performed in the generation of the basis vectors (equation (22a)) and the 
formation of the reduced matrices (equation (17)). Also, the back 

substitutions and formation of new reduced matrices required for augmenting 
the basis, or the back substitutions required to improve individual 
eigenmodes can represent significant contributions. 

4.4 Improvement of the Eigenvector Derivatives 

The eigenvalue and eigenvector derivatives are calculated from the 
reduced equations (24) to (29). Since the reduced basis is developed for 
representing the eigenvectors (and not their derivatives) and also since the 
reduced eigenvectors may not be accurate (the accuracy of the full 
eigenvector may have been improved), the eigenvector derivatives generally 
require some improvement. As discussed previously, one option would be to 
increase the characterization of the reduced basis. However, for the 


23 


eigenvector derivatives this is not considered viable for the following 
reasons: 

1. The reduced eigenvalues and eigenvectors have to be 
recalculated. 

2. In addition to reforming the reduced stiffness and mass 
matrices the derivatives of these matrices must also be 
calculated. 

3. Adding more vectors to the reduced basis increases the size 
of the reduced equations and the cost of obtaining their 
solution. 

4. The methods discussed previously for augmenting the reduced 
basis would not necessarily result in improved eigenvector 
derivatives. 

Two methods for improving the accuracy of the eigenvector derivatives are 
presented subsequently. 

4.4.1 Method Based on the Inverse Power Method 

The well known inverse power method for the generalized eigenproblem 
is given by the recursion relations 

([K] - V[M])(X.) <J + 1) . V)[H|{X.) <:|) 

X«> - 1 / ({X.] <j)T |M]{X 1 } <J) ) 1/2 

which results in the convergence of X^ to the eigenvalue closest in 

magnitude to the shift parameter \> and {X^} to the corresponding 

eigenvector. In the present study the eigenvalues and eigenvectors are 
already known to a required level of accuracy and iterations through 
equations (56) and (57) would result in very small changes. Noting this, 
equation (56) is differentiated with respect to an independent design 
variable and the superscripts on the eigenvalue, eigenvector and 

eigenvalue derivative (but not the derivative of the eigenvector) indicating 
the iteration are discarded. The resulting recursion relation is 


(56) 

(57) 


UK] - v[M] ) 


3{X.) 

~avT 


(j + 1) 


OX. 0[M] 

i[M] + X,- 


3V, 


i av, 


3[K] 

“avT 


(x.) 


9{X } 

+ (X.-v)[M]-^- 


(j) 


(58) 


If the full eigenvector was improved by iteration through equation (54) it 
is important that the derivative of the eigenvalue be calculated from the 
full equation (5) instead of the reduced equation (27) since the reduced 
eigenvector is less accurate. 

A major advantage for the procedure defined by equation (58) is that 
an existing decomposed coefficient matrix is used to avoid the cost of 
additional decompositions. This suggests an alternate development of 
equation (58). Returning to the full equation for calculating the 
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eigenvector derivatives (3a) it is noted that there is only a small 
difference between the matrix on the left-hand side ( [K] —X. [ M ] ) and the 

existing decomposed matrix ([K]-v[M]). The matrix on the left hand side can 
be split resulting in 


9{X.) 3{X } 

<[K] - V[M1)^- - 

k k 

{X.} (59) 

The second term on the left-hand side of (59) is then moved to the right- 
hand side. Using the current estimate for the eigenvector derivative, the 
right-hand side is evaluated and the equation solved for a new estimate of 
the eigenvector derivative. This relation is identical to that given by 
equation (58) and offers further insight into the procedure. 

4.4.2 Method Using Preconditioned Conjugate Gradient Iteration 

The preconditioned conjugate gradient (PCG) method for solving 
symmetric positive definite linear systems has been studied for many years 
[26,27]. In solving the system of equations 


U /V . ° l 1 1 J 

gy-[M] + X. 




3V, 


[A] {X} = {B} 


(60) 


the PCG method attempts to accelerate the convergence of the iterative 
process by choosing a preconditioning matrix [A] q such that 


[A] = I A ] o - [A] 


(61) 


The matrix [A] is the difference between the coefficient matrix [A] and the 
preconditioning matrix [A] q . At each iteration a system of equations with 

coefficient matrix [A] must be solved, so [A] should be chosen such that 

this system is easy to solve. In addition, the matrix [A] q should be chosen 

to in some way approximate [A]. In the present study, the splitting shown 
in (61) is motivated by the existence of a previously decomposed matrix. 

The objective of the PCG iteration is to solve equation (3a) for the 
derivative of the eigenvector using the approximate result obtained from 
equation (24) as an initial estimate. The selection of the preconditioning 
matrix ( [K] — v[ M ] ) as an existing decomposed matrix, results in the splitting 

UK]-X.[M]) = ( [K]-\>[M] ) - (X.-v) [M] 

The PCG algorithm as presented in [28] is as follows: 
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3. 

4. 

5. 

6 . 

7. 


3{X.} VU ' 

- ([K]-X.[M])-g^- , j = 0 (62) 

k 

Solve for the preconditioned residual {Y}^^: 

([K]-v[M]){Y} (j) = {R} (j) (63) 

Compute the orthogonalization coefficient, b., using 

T t 

bj = ([K]-v[M]){Y} (j) / {Y} (j ” 1) ([K]-v[M]){Y} (:i ~ 1) , j > 1, 

= 0* j = 0 (64) 

Update the conjugate search direction vector 

{Z} (j) = {Y} (j) + b.{Z} (j_1) (65) 

Compute the step length along the search direction, a., using 
X x 

a j = (Y} (j) ([K]-v[M]){Y} (j) / {Z} (j) ([K]-X.[M]){Z} (:i) (66) 

Update the solution 

3{X ) (j + 1) 3{X } (j) 

3V. = 3V, + a j* Z * < 67 ) 


Even though the PCG method is intended for solving positive definite systems 
and the present equations are singular, the nonsingularity of the 
preconditioning matrix allows the implementation of the procedure. The 
effectiveness of this algorithm is discussed in the succeeding section. 
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V. NUMERICAL STUDIES 


5.1 General 

This section describes the implementation of numerical algorithms 
for calculating free vibration sensitivity derivatives from reduced finite 
element equations. Their performance compared to the exact solution is 
described subsequently. Henceforth, the term "exact solution" is used to 
denote the solution obtained using the full system of equations of the 

structure. 

The structures considered in this study are composed of beam 

elements. The following three structures are considered: 

1. A non-symmetric triangular double-laced lattice beam ten 

bays in length, supported at one end (Figure (1)). 

2. An elastic line model of the Bell AH-1G Cobra helicopter 

(Figure (2)). 

3. An orthogonal tetrahedral lattice beam ten bays in length, 
supported at one end. The dimensions of each bay 
correspond to the current design of the keel beam of NASA's 
space station (Figure (3)). 

A simple program was written for modelling frame type structures 

using both displacement and mixed finite elements. The program was checked 

by comparing the free vibration results with those from the commercial 

finite element code EAL for the lattice beams and NASTRAN for the Cobra 

elastic line model. The primary motivation for developing this specialized 
code was to include the capability of assembling matrices of the derivatives 
of stiffness and mass coefficients of the structure and also to include 
mixed finite element models. Both of these tasks are not easy to implement 
in commercial finite element codes. 

The first structure considered, the triangular lattice beam, was 
developed as a test bed to provide a quick means of evaluating the 
algorithms presented. The cross-section is non-symmetric to further 
simplify the problem by eliminating repeated eigenmodes. The design 

variables are chosen to be the cross-sectional area (axial stiffness) and 
moment of inertia (bending and torsional stiffness) of the longerons in the 
last bay (next to the support, numbered 1, 2 and 3 in Figure(l)). 

The two additional cases represent more practical problems. The 
Cobra elastic line model is a simplified representation of the actual Cobra 
helicopter structure that is being used for sensitivity studies [28]. In 
addition to beam elements it contains linear spring and rigid body elements 
and uses a lumped approach in modelling the mass of the structure. The 
design variables are shown in Figure (2) and include the stiffness of linear 
spring and beam elements and certain concentrated masses. The other case 
considered is a ten bay section extracted from the support structure of the 
current proposed configuration of the NASA space station. The design 

variables are very similar to those chosen for the triangular lattice. 

5.2 Comparisons of the Exact and Approximate Vectors 

In comparing vectors obtained from the reduction method 
(eigenvectors and their derivatives) to those obtained from the solution of 
the full equations the following error norm is used [14]: 
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N 

e = - 2 [ (X -X ) / 

N i=l 11 


|X | ]‘ 
1 max 1 1 


( 68 ) 


{X} is the exact vector and {X} is the vector obtained from the reduction 

method. IX I is the maximum absolute value of the exact vector and N is 
max 

the number of components. This norm is similar to the r.m.s. error norm but 
is weighted to give less importance to the smaller components of {X} which 
are usually of less practical importance. 

It is important to distinguish between the norm presented here using 
the exact solution and the one presented in the previous section using a 
residual from the full equations. Experience in this study has shown that 
values for the exact error norm of 0.001 or less indicates an accurate 
approximate vector. This corresponds to an error for the largest component 
of a few percent or less. Values for the residual error norm of 0.1 or less 
show a similar accuracy. On the other hand, values for the exact error norm 
of 0.01 or greater indicate that the approximate vector shows very little 
resemblance to the exact. For the residual error norm values of 1.0 or 
greater indicate poor accuracy. 

In addition to the error norms, portions of the exact and 
approximate vectors are plotted. The components and nodes for which the 
vectors are plotted are chosen by examining the eigenvectors for the largest 
displacement and rotation. These plots provide a better physical indication 
of the accuracy of the vectors. 

5.3 Verification of the Results 


To insure the reliability of the presented results they are checked 
by independent programs using finite difference operators to evaluate the 
derivatives. For the full equations, eigensolutions are calculated for a 
positive and negative perturbation of each design variable from the 
equations 

(|K]< ViflV >- (X.)< ViSV >= 0 (69) 

where the superscript (V±AV) indicates the positive and negative 
perturbation respectively. The eigenvalue and eigenvector derivatives are 
then calculated from the first order central difference formulae: 


ax. x< v+ ®>- x‘ v - lv > 

1 _ 1 1 

3V ~ 2AV 

3{X.} (X.}< v+Av >- {X.} (V - AV > 

“av - ~ 2AV 


(70) 

(71) 


In the development of the reduced equations it is assumed that the basis 
vectors are invariant with respect to the design variables. Thus in the 
finite difference approximations the basis vectors are determined for the 
unperturbed structure and used to reduce the perturbed equations resulting 
in 
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[K1< V44V) - [ri< v ) T iK]< Vi6V >(r]< v > 

[R] (v±iv) . m <v)T [H]< v± 4 V >[r] <v) 

(IKl< Viiv >- X< ViAV >[Rl <Vi4V >) (*.]< Vi4v > 


The approximate derivatives are then calculated from the relations 


ax. 

i 

^(V+AV) 

i 

av ~ 

2AV 

3{X.} 

= m <v > 

av 




(V+AV) 


c* ± > 


(V-AV) 


2AV 


(72a) 

(72b) 

(73) 


(74) 

(75) 


In all cases very close agreement was found between the analytic and finite 
difference results. 


5.4 Accuracy of the Derivatives Predicted by the Reduction Methods 

In the present study both static condensation and generalized 
dynamic reduction are investigated as a means to reduce the governing 
equations for the displacement model of the triangular lattice beam. The 
results for the lowest eight vibration modes are presented in Tables (1) and 
(2) and the first four eigenvectors and their derivatives are plotted in 
Figure (4). 

In the application of static condensation to the triangular lattice 
beam, 36 out of the total 180 degrees of freedom are retained in the 
analysis. These correspond to the six degrees of freedom at each of the 
three nodes on the first and third cross-sections of the beam, numbered 
sequentially from the free end. Only the lowest two eigenvalues and their 
corresponding eigenvectors are reasonably accurate. For higher modes the 
accuracy of the eigensolution falls off very rapidly. The derivatives of 
the eigenvalues are less accurate than the eigenvalues themselves, but are 
within ten percent of the exact value. Unfortunately, the derivatives of 
the eigenvectors show very little or no resemblance to the exact solution, 
even for the lowest two modes. 

With generalized dynamic reduction a much higher degree of reduction 
is possible. Several choices for the number of basis vectors and iterations 
were tried and the results indicated that accurate modes could be obtained 
with as few as six basis vectors. Attempts at using a larger number of 
vectors resulted in slightly more accurate eigensolutions but had little 
effect on the accuracy of the eigenvector derivatives. The cost of reducing 
and solving the reduced equations increases rapidly with the number of basis 
vectors so it is generally more efficient to use the smallest number of 
basis vectors possible. Although fewer modes are obtained from the reduced 
set of equations and more shifts are required, it is still generally more 
efficient . 

The results shown in Figure (4) and Tables (1) and (2) are obtained 
by using six basis vectors that are generated with four iterations on the 
initial random vectors. In each reduction the lowest four of the calculated 
modes are retained. The shift value is initially chosen to be zero and 
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calculated from the estimates of the higher eigenvalues on subsequent 
reductions. Although only the first eight modes are shown, the first twelve 
were calculated accurately. Even with such a high degree of reduction, the 
eigensolution is obtained much more accurately than with static 
condensation. The accuracy of the eigenvalue derivatives is also 
considerably better although they are still less accurate than the 
eigenvalues. The eigenvector derivatives in some cases resemble the exact 
results but for the most part, the reduced basis does not adequately 
characterize the eigenvector derivatives. 

5.5 Implementation of the Mixed Method 

In the previous section the triangular lattice is modelled with a 
displacement beam element. This section investigates the modelling of an 
identical structure with a mixed beam element. The element has 3 
displacement nodes and 2 stress nodes resulting in quadratic approximation 
of the displacement field and linear approximation of the stresses. In 
contradistinction, the displacement element has two nodes and uses 
interpolation functions developed from classical beam theory that result in 
quadratic approximations for the axial displacements and cubic 
approximations for the bending displacements [30]. Also, the mixed element 
includes shear deformation and the displacement element does not. The 
results for the lowest four modes of the mixed model are presented in Table 
(3) and Figure (5). Comparing these with those obtained for the 
displacement model it is seen that the full equations yield almost identical 
results . 

The full mixed model contains 900 displacement and 240 stress 
degrees of freedom. Through the use of generalized dynamic reduction the 
total number of degrees of freedom is reduced to just 6. Again comparing 
the results from the two formulations it is seen that the reduced mixed 
model does not provide quite as accurate results as the reduced displacement 
model of the previous section (reduced from 180 to 6 degrees of freedom) but 
this is probably due to the large difference in the number of degrees of 
freedom. However, the results are very similar in that the eigensolution 
and the eigenvalue derivatives are determined fairly accurately but the 
eigenvector derivatives are not. It is apparent from both the displacement 
and mixed models that some type of improvement is necessary for the 
eigenvector derivatives obtained from the reduction method. 

The advantages of this method result from the simplified form of the 
governing equations. As discussed in a preceding section, this simplifies 
the calculation of the derivatives of the stiffness matrix when the design 
variables are the inverse of a cross-sectional area or moment of inertia. 
Also, shear deformation can be included in the formulation very easily. 
This study shows that the mixed model is capable of providing the same 
results as the displacement model. 

5.6 Iterative Refinement of Eigenvector Derivatives 

Two algorithms are implemented for improving the accuracy of the 
eigenvector derivatives that were obtained using generalized dynamic 
reduction on the displacement model of the triangular lattice beam. These 
algorithms are presented in the preceding section. The first method is 
based on the inverse power method and will be referred to as IPM while the 
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second uses a preconditioned conjugate gradient algorithm and will be 
referred to as PCG. 

In the present implementation the IPM or PCG iterations can be 
terminated by three separate criteria. A maximum of two iterations is 
allowed, but if the residual error norm for the vector drops below 0.1 or is 
greater than on the previous iteration the procedure terminates. The 
results of using IPM and PCG iteration are shown in Tables (4) and (5) 
respectively for the lowest eight modes. For both IPM and PCG, the results 
are plotted in Figure (6). It is seen from these results that one iteration 
with either method in most cases provides a dramatic improvement in the 
accuracy of the approximate vector. In some cases a second iteration is 
required to meet the convergence tolerance but in a few of these cases the 
second iteration diverges (Table (5), PCG iteration, derivative of 
eigenvectors 6-8 with respect to also, not shown, IPM iteration, 

derivative of eigenvector 12 with respect to l^). This may be a result of 

the near singularity of the matrix on the left-hand side of the governing 
equations. Further investigation with more iterations shows an oscillatory 
behavior for all the derivatives but this does not usually occur in the 
first two iterations. It is observed to occur after fewer iterations for 
cases in which the approximate eigensolution is not quite as accurate, but 
occurs for all modes after several iterations. Although the cases for which 
oscillation occurs on the second iteration do not quite meet the convergence 
tolerance, examination of the plots in Figure (6) show that the approximate 
vectors are reasonably accurate. 

From the examination of Tables (4) and (5), it appears that the IPM 
iterations are slightly more effective at improving the accuracy of the 
eigenvector derivatives than PCG, since it generally leads to smaller values 
of the error norm. However, the plots of the vector shown in Figure (6) 
show the difference to be insignificant. In terms of efficiency, PCG has an 
advantage since the residual needed for the convergence check is a by- 
product, whereas in IPM iteration it must be calculated separately. 

A topic of primary interest is the efficiency of the reduction 
method as compared to the solution of the full system of equations. The 
programs used in the current implementation are not optimized sufficiently 
to make a comparison of execution times meaningful. Also, the structures 
investigated are smaller than the intended application due to the central 
storage limitation (130560 words) of the CYBER 175 computer used in the 
analysis. In order to assess the efficiency, an accounting of operations is 
performed. The comparison is made between the solution of the full system 
of equations by Nelson's method (equations (5)-(7)) and the solution of the 
reduced equations with subsequent iterative refinement (1 iteration) by the 
PCG method. With two major exceptions, the operations required in both 
cases are almost identical. The exceptions are the decomposition of a 
coefficient matrix for each eigenvalue in the solution of the full equations 
and the reduction of the derivatives of the stiffness and mass matrices to 
obtain the reduced equations. Both of these have no counterpart in the 
other method. Considering the usual sparsity of the derivatives of the mass 
and stiffness matrices, the latter is usually much less time consuming. 
Since the decompositions usually account for a large percentage of the 
computer time the reduction method can result in significant savings. 
However, if the design variables affect large portions of the structure or 
there are an unusually large number of design variables, the solution of the 
full equations may be more attractive. 
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5.7 AH-1G Cobra Elastic Line Model 


The Cobra elastic line model (developed by Hanson and Mur thy [29]) 
is a very coarse representation of the actual helicopter structure as shown 
in Figure (2). The displacement model contains 336 degrees of freedom 
before the inclusion of rigid body elements and the application of boundary 
conditions. The rigid body elements reduce the number to 252 and the 
boundary conditions (chosen to eliminate the rigid body modes) further 
reduce the number to 246. 

The lowest eight vibration modes (shown in Figure (7)) can be 
divided into groups according to the portion of the helicopter to which they 
correspond. Modes 1 and 2 represent rocking of the main rotor pylon 
assembly; 3, 4, 7 and 8 correspond to modes of the aft fuselage, tailboom 
and vertical fin; and modes 5 and 6 to the front fuselage section. Two 
design variables are defined within each of these sections and are marked by 
circled number in Figure (2). For the main rotor pylon the design variables 
are the stiffness of the springs supporting the pylon in the vertical 
direction^ (T) ) and the bending stiffness of the upper section of the rotor 
mast ( (2) ). For the aft fuselage modesthe bending stiffness of a small 
section is chosen as one design variable ((3)) and the mass of the tail 
rotor as the other ( (4) ). Likewise for the front fuselage the bending 
stiffness of one element is chosen as a design variable ( (|) ) and also the 
mass of the gun turret ((6)). 

The results obtained from the full equations (exact) and using 
generalized dynamic reduction with 12 basis vectors are presented in Figure 
(7) and Table (6). The basis vectors are obtained by performing 4 
iterations on 6 random and 6 unit vectors. The unit vectors contain a 
single unit component in one of the locations corresponding to the masses of 
the main and tail rotors. This ensures that these masses are excited and 
improves the accuracy of the associated vibration modes. The results show 
that this provides good accuracy for the eigensolutions. 

In examining the accuracy of the eigenvalue derivatives some very 
large errors are discovered, however these correspond to modes that are 
effectively invariant with respect to the particular design variable. Since 
the reduction method does find the eigenvalue derivatives to be relatively 
small, these large errors do not pose a problem. For efficiency logic 
should be included in the program to avoid computation of the associated 
eigenvector derivatives. 

The eigenvalue derivatives of importance are predicted accurately. 
In addition, a few of the eigenvector derivatives are also reasonably 
accurate, but the majority are not. The application of PCG iteration proved 
very effective at obtaining accurate eigenvector derivatives. 

5.8 Orthogonal Tetrahedral Lattice Beam 

In addition to the previously presented triangular lattice that 
represents a specialized problem, it is desired to present a repetitive 
beam-like lattice structure with some practical significance. The keel beam 
of the proposed NASA space station consists of several orthogonal 
tetrahedral lattice beams (see Figure (3)) and for this study a ten bay 
length is investigated. 

The beam is modelled with a displacement element and contains 240 
degrees of freedom. The model is reduced using generalized dynamic 
reduction with eight basis vectors. For the initial reduction (shift = 0.0) 


32 



the basis vectors are obtained by performing 4 iterations on eight random 
vectors and the lowest two eigenmodes are retained. On subsequent 
reductions 6 iterations are performed on the random vectors and two modes 
are retained. The results for the first eight modes are presented in Figure 
(8) and Table (7). It is observed that for both the present case and the 
triangular lattice, the minimum number of basis vector required is l/30th 
the total , number of degrees of freedom, however the present case requires 
more iterations and yields a smaller number of accurate modes per reduction. 
The first twelve modes were calculated accurately. It is suspected that 
more could be calculated, however this was not attempted. 

The derivatives of the eigenvalues are obtained accurately from the 
reduction method and PCG iteration is used to improve the eigenvector 
derivatives. Using the same termination criteria discussed previously, no 
cases are present in which the solution diverges in the first two 
iterations. The PCG method again demonstrates good effectiveness at 
improving the accuracy and in most cases only one iteration is required. 
However, comparing the error norms in Table (7) with the plots in Figure 
(8), one peculiarity is discovered. For the derivative of eigenvector 3 
with respect to the exact error norm is 0.0016 after PCG iteration which 

indicates a reasonably accurate vector. The plot of the vector shows good 
agreement for the rotations but poor agreement for the displacements. This 
is most likely due to the difference in ,the size of the components 
(approximately two orders of magnitude) and brings out a small weakness of 
the procedure. In examining the remainder of the results a few more cases 
are discovered in which the smaller terms are not calculated accurately. 
However, in many practical applications these components are of little 
importance. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


6.1 Summary 

In the present study computational procedures are presented for 
calculating the sensitivity derivatives of the vibration frequencies and 
eigenmodes of framed structures. The procedures incorporate two key 
elements. The first is the use of dynamic reduction techniques to 
substantially reduce the number of degrees of freedom. The reduction 
techniques considered are static condensation (Guyan's reduction) and a 
generalized dynamic reduction technique. The second key element is the 
application of iterative techniques to improve the accuracy of the 
derivatives of the eigenmodes. The two procedures presented are based on 
the inverse power method and the preconditioned conjugate gradient 

technique. 

Numerical studies are performed on three frame type structures using 
displacement and mixed finite element models. Error norms are introduced to 
assess the accuracy of the eigenvalue and eigenvector derivatives obtained 
by the reduction techniques. 

6.2 Conclusions 

On the basis of the numerical studies performed in this report the 
following conclusions seem to be justified: 

1) In the structures considered, generalized dynamic reduction allowed 

the reduction of the equations to a small number of degrees of 
freedom and provided a decomposed coefficient matrix useful in 

improving the accuracy of the eigenvector derivatives. Static 
condensation was sensitive to the choice of the free degrees of 

freedom. Neither method was capable of reliably predicting the 

derivatives of the eigenvectors. 

2) The mixed and displacement models used in the present study resulted 
in eigensolutions and their derivatives that were in very close 
agreement. Also, similar accuracy was obtained from the reduced 
models in both formulations. The advantage of the mixed formulation 
results from the simplified calculation of the derivatives of the 
stiffness matrix with respect to design variables that are the 
inverse of sizing quantities. Also, the simple form of the 
governing equations allows shear deformation to be included easily. 

3) Both the exact and residual vector error norms provide a good 

indication of the accuracy of the approximate vectors. A good 
correlation was observed between the two norms, the error of the 

largest component and the plots of the vectors. 

4) Both the PCG and IPM procedures for improving the accuracy of the 
approximate eigenvector derivatives proved very effective. In most 
cases one iteration was sufficient to obtain very accurate results. 
The efficiency of these procedures compared with other solution 
methods using the full equations depends on several factors. For 
cases in which the decomposition of the coefficient matrices results 
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in a large percentage of the full system solution time and the 
derivatives of the stiffness and mass matrices are sparsely 
populated, significant savings can be realized. 

6.3 Recommendations 

The following recommendations are presented in the hope that they 

will provide direction for future research in the topics investigated in 

this study: 

1) The effectiveness of the proposed procedures should be investigated 
for structures containing repeated or very closely spaced 
eigenvalues. One possible way of doing this is discussed in section 
two but numerical studies are not performed. 

2) The procedures presented in this report should be implemented in a 
commercial finite element code and numerical studies performed on 
several very large and complex structures. The size of the 
structures investigated in this study was limited by the central 
storage capacity of 130,560 words for the CYBER 175 computer. 

3) Detailed studies of the numerical behavior of the IPM and PCG 
procedures should be performed. The methods showed a tendency to 
oscillate after several iterations and the causes of this should be 
identified. 

4) The effectiveness of using static condensation to provide starting 
iteration vectors for the generalized dynamic reduction should be 
investigated. 
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Inertia, Axis 1 
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Inertia, Axis 2 

Torsion 
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Longeron 

A 

x i 

X 2 

J 

Batten 

0.5A 

0.1083I 1 

0.1083I 2 

0.1083J 

Diagonal 

0.5A 

0sl083I 1 

0.1083I 2 

0.1083J 


A - 3.0 X 10 -5 m 2 E = 6.895 X 10 10 N/m 2 

Ij - I 2 - 6.0 X 10" 9 m 4 G “ 2.652 X 10 10 N/m 2 

J - 1.2 X 10 -8 m 4 p - 2768 kg/m 3 

V 2 » Cross-sectional area of member s(^, (5) and(3) 

Vj = Moment of inertia and torsion constant of member^^ , ^^and(^ 
c) Design Variables and Material Properties 

Fignre (1) Dn symmetric triangular double laced beamlike lattice used in the 
present study. 
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a) Actual helicopter airframe structure. 



b) 


Vj = stiffness of linear elastic elements © ■ 7.881X10"’ N/m 
V 2 = moment of inertia of beam element © = 4.995X10 -5 m 4 


Vj = moment of inertia of beam element © ■ 4.162X10 - ^ m 4 
Ij = 1.167 V 3 , I 2 = 1.39 V 3 

V 4 * rigid mass © “ 13.97 kg 

Vj = moment of inertia of beam element (J) *■ 4.162X10 m 4 

Ij * 1.471 V $ , I 2 * 1.275 V $ 

* 

V 6 = rigid mass (b) - 114.9 kg 

E - 6.895X10 9 N/m 2 G - 2.652X10 9 N/m 2 

Design variables and material properties. 



Figure (2) Elastic line model of the Bell AH-1G Cobra used in the present 
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Table 1. Accuracy of the eigenvalues and their derivatives. 
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T.bl« 1. and 2. Reaulta for the displacement model of the ten bey triangular lattice bean (Figure 1.). The approximate re.ulte are obtained using atstic 
condensation with 36 retained dof, and also generalized dynamic reduction with 6 basis vectors. 
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Note: 

Table 4. and 5. Results of iterative refinement for the displacement model of the ten bay triangular lattice beam (Figure 1.). The 
approximate results are obtained using generalized dynamic reduction with 6 basis vectors. The derivatives of the 
eigenvectors are refined by performing the indicated number of iterations. The exact and residual error norms are 
calculated as shown in equations (68) and, (52) and (53) respectively. 
















Table 6a. Accuracy of eigenvalues and eigenvectors 
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Tables 6. Continued on next page 























Derivative of Eigenvalue and Derivative of Eigenvalue and 

Eigenvector w.r.t. V- Eigenvector w.r.t. 



78 


Note; 

Tables 6. Results for the displacement model of the Cobra elastic line model (Figure 2.). The approximate results are obtained using 
generalized dynamic reduction with 12 basis vectors. The derivatives of the eigenvectors are refined using preconditioned 
conjugate gradient iterations. The exact error norm is given by equation (68) and compares the exact and approximate vectors. 













Table 7a. Accuracy of eigenvalues and eigenvectors. 
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Note : 

Tables 7. Results for the displacement model of the orthogonal tetrahedral lattice beam (Figure 3.). The approximate results are obtained 
using generalized dynamic reduction with 8 basis vectors. The derivatives of the eigenvectors are refined using preconditioned 
conjugate gradient iterations. The exact error norm is given by equation (68) and compares the exact and approximate vector. 
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